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We develop a semi-classical field method for the study of the weakly interacting Bose gas at finite 
temperature, which, contrarily to the usual classical field model, does not suffer from an ultraviolet 
cut-off dependence. We apply the method to the study of thermal vortices in spatially homogeneous, 
two-dimensional systems. We present numerical results for the vortex density and the vortex pair 
distribution function. Insight in the physics of the system is obtained by comparing the numerical 
results with the predictions of simple analytical models. In particular, we calculate the activation 
energy required to form a vortex pair at low temperature. 
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I. INTRODUCTION 

Classical field theories are a widespread and flexible 
tool for the study of many aspects of the physics of ultra- 
cold Bose gases. Developed in particular to address time- 
dependent problems related to dynamical aspects of the 
Bose-Einstein phase transition [JlH they can also be used 
to study thermal equilibrium properties of the weakly in- 
teracting Bose gas [1, 0] . A major example is the quan- 
titative prediction of the shift of the Bose-Einstein con- 
densation temperature due to atomic interactions [5( that 
has been obtained by means of a Monte Carlo sampling 
of a classical field model [3| : as long as the physics of the 
system is determined by the low-energy modes, classi- 
cal field models provide reliable results on the full quan- 
tum problem. Classical field techniques have also been 
applied to obtain analytical and numerical predictions 
for reduced dimensionality Bose systems 0, 0, H, Q , in- 
cluding the calculation of the critical temperature for the 
Berezinskii-Kostcrlitz-Thouless transition in two dimen- 
sions [ToL Hlj |. However, an ultraviolet cut-off has to be 
introduced in most of these classical field techniques in 
order to avoid ultraviolet divergences analogous to the 
blackbody catastrophe of classical statistical mechanics, 
and this raises the problem of a possible cut-off depen- 
dence of some of the physical results. 

On the other side, several exact reformulations of the 
many boson problem have been developed. Although 
they have successfully served as a starting point for 
Quantum Monte Carlo simulations [H, EH of the ther- 
mal properties of Bose systems such as liquid Helium and 
ultracold atomic gases 0, EE E3] , they often lack the in- 
tuitiveness of classical field theories where the physics is 
described in terms of a simple distribution function in 
the functional space of c-number fields. 

The present paper is devoted to the development, the 
validation, and the first application of a semi-classical 
field theory which tries to combine a regular behavior in 
the ultraviolet limit with a transparent intuition of the 



physics of the system. As in classical field theories, the 
density matrix of the Bose system is written in terms of a 
distribution in the space of c-number fields. In the semi- 
classical theory, this distribution is however much more 
complex than a simple Boltzmann factor exp(— E/ksT), 
where E would be the Gross-Pitaevskii energy of the field 
configuration, and has to be obtained as the result of an 
imaginary-time Gross-Pitaevskii evolution starting from 
an initially uniform distribution in functional space. 

A first application of the method is then presented to 
the study of thermal vortices in a homogeneous two- 
dimensional Bose gas, in particular their density and 
their pair distribution function. Experimentally, the 
two-dimensional Bose gas has been realized some time 
ago [U EH, but it is only recently that several exper- 
iments have given indications of the presence of vor- 
tices in finite temperature samples [13, 0, [2l[ , and this 
raises the question of the link between observable quan- 
tities (e.g the vortex density), and theoretical concepts 
such as the Berezinskii-Kosterlitz-Thouless (BKT) tran- 
sition [13, [U miH, Hi| . Most of the existing theoretical 
treatments neglect all density fluctuations other than the 
ones in the vicinity of a vortex core, and eventually map 
the 2D Bose gas problem onto the XY model of statistical 
mechanics [25| ■ Although this approximation is expected 
to provide a good description of atomic gases trapped in 
2D optical lattices [111, 0, [27J , it seems far from being 
accurate for spatially continuous systems: at tempera- 
tures of the order of the BKT transition temperature, 
the amplitude of the density fluctuations in the gas is 
not negligible as compared to the density itself [28] . Our 
work aims at going beyond this approximation so to fully 
include the effect of density fluctuations. The fact that 
it is based on c-number fields gives to the present semi- 
classical method an advantage over standard Quantum 
Monte Carlo techniques in view of the study of vortices. 

The paper is divided in two main parts. In the first 
part (Sec|TT|, we introduce the semi-classical method 
in the grand-canonical (Sec lII A| ) and in the canonical 
fSec lII Cl) ensembles, and we characterize its range of 
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applicability (Sec lII B| ). In the second part (Sec lIII[) . 
we discuss the physics of the two-dimensional Bose gas. 
The numerical results are presented in Sec lIII Al differ- 
ent observables are considered, e.g. the normal and non- 
condensed fractions, the density fluctuations, the vortex 
density, and the vortex pair-distribution function. In 
Sec lIII Bl the effect of Bose condensation on the vortex 
density in the finite size ideal gas is discussed analyt- 
ically; this requires the use of the canonical ensemble, 
which introduces new features with respect to the well- 
studied grand canonical case [29, 30]. In Sec lIII Cl a sim- 
ple model including the interacting case is developed to 
understand the numerical results, principally the ones for 
the vortex density an activation law of the form 

n Vl + oc exp(— A/fceT) is found in the low-temperature 
regime, and the dependence of A(T) on the system pa- 
rameters such as the interaction strength and the system 
size is discussed: the main qualitative differences between 
the ideal and the interacting gas behaviors are pointed 
out. Conclusions are finally drawn in Sec lIVI 



II. THE SEMI-CLASSICAL METHOD 
A. In the grand-canonical ensemble 

Consider a Bose field defined on an square lattice of 
TV points with periodic boundary conditions; V is the 
total volume of the quantization box and dV = V/ftf 
is the volume of the unit cell of the lattice. The Bose 
field operators ^(r) obey the Bose commutation relations 
$(r),&(r')]=5 r , r ,/dV. 

The state of the Bose field is described by the den- 
sity operator p, which can be expanded in the so-called 
Glauber-P representation on coherent states: 



p = VipP[ip)\coh:if>)(coh:i(>\, 



(1) 



where P[if>], the Glauber-P distribution, is guaranteed to 
exist in the sense of distributions but in general is not a 
positive nor even a regular function [3l|, HH, [H[ . tjj(r) is 
here a c-number field defined on the lattice, the coherent 
state is defined as usual as: 



|coh : ip) = exp 



-~M\ 



e X p|^dF^(r)#t( r )j | ), 

where 1 | 2 = dVJ2 r \' l l J { r )\ 2 > anc ^ * ne functional integra- 
tion is performed over the value of the complex field at 
each of the Af sites of the lattice: 

VtJj = Y[ dRe ty>(r)] dim [ip(r)] . (3) 

r 

The homogeneous Bose gas is described by the follow- 



ing second-quantized Hamiltonian: 

~h 2 k 2 
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^^rfF*t(r)^(r)*(r)*(r). (4) 



The single-particle dispersion relation within the first 
Brillouin zone is taken as parabolic with mass m, p is 
the chemical potential, and the interactions are modeled 
by a two-body discrete delta potential of strength go- 

The gas is assumed to be at thermal equilibrium at a 
temperature T, so that the unnormalized density opera- 
tor is p cq (/3) = cxp[— f3H] with f3 = l/fc^T. This density 
operator can be obtained as the result of an imaginary- 
time evolution: 



dp, 
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Pc q H] (5) 

(3, starting from the 



during the "time" interval r = 
identity operator p e q(r = 0) = 1. 

In the Glauber-P representation, the imaginary-time 
evolution takes the form of a Fokker-Planck-like partial 
differential equation: 
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(6) 



for the distribution function P[ip] in the phase-space of 
the c-number fields defined on the lattice. The derivatives 
with respect to the complex field ip(r) arc defined as usual 
as: 



d, 



(7) 



The first term in the right-hand side of (O acts on the 
weight of the wavefunction ip and involves the mean-field 
energy of the complex field ^>(r): 



em = J2 dv ^*( r ) fa - d + f E dv i^( r 



(8) 

ho is a shorthand for the single-particle Hamiltonian, 
whose fc-space form is ho — h 2 k 2 /(2m). 

The second term is a drift term consisting of the 
imaginary-time Gross-Pitaevskii evolution: 



mi*) 



i 



2dV 



1 



[ho-v + goH(r)\ 2 ] ^(r). 

(9) 



Finally, the diffusion terms involving the second-order 
derivatives are local in space, but have a non-positive- 
definite diffusion matrix: 



D(r) 



go 

4dV 



^ 2 (r) 
ip* 2 (r) 



(10) 
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A complete solution of the partial differential equation 
would provide the exact result of the lattice quantum 
field problem defined by the Hamiltonian (j4|). Unfortu- 
nately, the non-positive-definite nature of the diffusion 
matrix (fT0|) prevents the Fokker-Planck-like equation ([6]) 
from being directly mappable on a stochastic field prob- 
lem for tp. Some approximation schemes are therefore re- 
quired in order to perform numerical simulations within 
the Glauber-P framework. 

In our previous work Q , the high-temperature physics 
of the one-dimensional Bose gas was studied by keep- 
ing only the first term in the right-hand side of j6|. 
The resulting distribution in the phase-space of the 
c- number fields is the usual Boltzmann one P[ip] = 
exp(— E[ip]/kBT) in terms of the mean-field energy {8j. 
A better approximation is obtained by keeping also the 
drift force ([9]) and neglecting the diffusion term (fl0|) only. 
In this case, the partial differential equation ([5]) can be 
mapped onto a deterministic evolution for the field ip( r ) 
and a weight W: 



d T ip(r,T) 
d T W{r) 



1 



go |iMr,r)| 2 ] V(r, ' 
E[4>(t)]W(t). 



(11) 
(12) 



Physical quantities are then obtained as averages over 
the initial values for ip. A possible representation of the 
initial state p cq (r = 0) = 1 is to take the initial value 
of the field ip(r,T = 0) at each lattice point as uniformly 
distributed in the complex space and to take a constant 
initial weight W(r = 0) = wq. This leads to the semi- 
classical approximation for the density operator at tem- 
perature T: 



solved and their prediction compared to the exact quan- 
tum results. 

By defining the operators Ck,+ = (&k + Q-k)/v / 2 and 
c\ = (cik — fl-k)/(*V / 2), the Bogoliubov Hamiltonian 
(|14[) can be rewritten as a sum of terms involving inde- 
pendent k modes: 



Ti-Bo 
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CwCk.e 



(15) 



In this way, the Glauber-P distribution factorises as a 
product of independent factors involving the different k 
modes. To avoid double-counting of modes, the primed 
sum is restricted to those k vectors which are contained 
in an (arbitrarily chosen) half-space. 

Each term of the Hamiltonian (|15|) has the simple 
structure of a one-mode squeezing Hamiltonian: 



(16) 



with the kinetic energy coefficient E k = f?k 2 /(2m) and 
the c operator corresponding to any of Ck,± in (fTS"]) . 

Since the Hamiltonian (|16p is quadratic, the exact 
Glaubcr-P distribution for the thermal equilibrium state 
can be analytically obtained by means of standard tech- 
niques [33j |. as well as its semi-classical approximation: 
as shown in the Appendix [XJ both distributions have a 
Gaussian form, 



P( 7 )cxe-( Rc ^ 2 / CT V( Im ^ 2 /^ 



(17) 



Psc = J Zty(0)W()8)|coh : VX/3))(coh : V(/3)|, (13) 

where both W(/3) and ip((3) depend on the initial value 
of the field V>(0). 

As the diffusion term (flU)) is proportional to the in- 
teraction strength go, the semi-classical approximation 
becomes exact in the case of the free Bose field, i.e. for 
an ideal Bose gas. As a consequence, it does not suffer 
from the typical ultraviolet divergences of classical field 
theories, even in presence of interactions. 



B. Limits of validity 

In order to validate the semi-classical approximation 
and appreciate its power and its limits, it is interesting to 
apply it to the simple case of the Bogoliubov Hamiltonian 
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This Hamiltonian being quadratic in the field operators, 
the semi-classical equations (|1HI12[) can be analytically 



The widths for the exact distribution are given by 



("?). 




cotanh ( — — | — 1 

2 ; 



cotanh 



(3e k 



1 



(18) 
(19) 



where e k = [Ek(2/i + E^)] ' is the energy of the Bogoli- 
ubov mode. When the temperature is too low, (cr^) ex be- 
comes negative, so that the Glauber-P distribution ceases 
to exist as a regular function [32, HH . The correspond- 
ing lower bound on the temperature is plotted in FigfT] 
Two limiting cases are easily isolated: for low-energy 
modes such that E k — » 0, the positivity condition for 
the Glauber-P distribution is the simple one fc^T > /i. 
For high energy modes, the condition is instead more 
stringent, k B T > (E k + fi)/\ log(fi/2E k )\. 

The widths for the semi-classical approximation are 
given by 
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,0(E h +2,j,) 
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FIG. 1: In the Bogoliubov model, minimal value of the 
temperature T m i n ensuring regularity and positivity of the 
Glauber-P distribution in a mode k, as a function of the ki- 
netic energy coefficient Ek of the mode. 
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FIG. 2: (Color online) In the Bogoliubov model, mean energy 
in a mode as a function of the mode kinetic energy coefficient 
Ek for different values of the temperature ksT / u = 0, 2, 3, 5 
(from bottom to top). Solid lines: quantum result. Dashed 
lines: semi-classical theory. Dotted lines: classical field ap- 
proximation. 



As expected, they remain positive at all temperature. 

These results are the starting point for detailed com- 
parison of the semi-classical predictions to the exact 
quantum results for the most significant observables. Let 
us start with the mean energy. The semi-classical value 
is: 



(H 



1 sc 



E k 



Ek + 2/i 
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which is to be compared to the exact value 
ek . e/c - {Ek + m) 
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af3e k 
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(22) 



(23) 



An order by order comparison can be performed in the 
high-temperature limit by expanding (|22|) and (|23[) in 



powers of (3: 
(Hi)sc 

(Wl)ex 



kl;T _^k+Ji +0 ^ {Ek+2 ^ (24) 
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E k + t i 



0(l3e 2 k ). 



(25) 



Agreement is found not only on the classical term k B T, 
but also on the subleading constant term — (Ek + m)/2, 
which would instead be missed by a simple classical field 
theory. 

A more detailed comparison is obtained by working 
out two limiting regions. In the low energy limit, one has 



lim (Hi) ex. = k B T - — 



(26) 



lim (Hi)sc = k B T-'^+ -8^ + 0(/?V) : (27) 

the relative error of the semi-classical result is therefore of 
the order of (/3/x) 2 /6, i.e. very small provided fcgT ;g> ii. 
In the high energy limit e k — > oo, one has instead [3J| 



(H 



l/SC 



cosh(/3/Lt) e k e-^ k . 



(28) 



In the high temperature regime where cosh(/3/i) ~ 1, this 
semi-classical prediction almost coincides with the exact 
value (f2"3"]) once the zero-point energy is subtracted from 
the quantum value. This shows that the semi-classical 
theory does not suffer from any ultraviolet divergence 
coming from the zero-point energy, nor from the typical 
black-body catastrophe of classical field theories. 

In summary, the semi-classical theory is able to accu- 
rately reproduce the value of the average energy under 
the assumption that the temperature is higher than the 
chemical potential, k B T ^S> /i. Examples of plots of the 
mean energy of the different Bogoliubov modes as a func- 
tion of Ek are presented in Figj2] for the semi-classical 
theory, the classical field approximation, and the exact 
result. The agreement of the semi-classical theory with 
the exact result is already remarkable for temperatures 
only a few times higher than the chemical potential, while 
the classical field approximation is quite crude in predict- 
ing a constant mean energy k B T independent of the mode 
energy. 

Another observable that we consider is the normal 
fraction f n , defined as 



fn 



(Pi 



Nmk R T 



(29) 



where P x is the x component of the total momentum of 
the system. This quantity f n estimates the response of 
the Bose system to a gauge field, e.g. a magnetic field in 
the case of c harg ed particles, or a rotation in the case of 
neutral ones [35|, [36J . 

The exact quantum result of the Bogoliubov theory 
1371 has the form 



(Pi 



^2h 2 k 2 x n k (n k + l) 



(30) 
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where = {e^ tk — is the quantum mean occupa- tribution function, 
tion number of the Bogoliubov mode. The semi-classical 
approximation is instead given by 



<? ( V 



|rt( r ) *t( r ')^( r ')^( r ) 



(33) 
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sc 



'sc 



(31) 



It is interesting to compare the expression between 
square brackets to the quantum value nk{nk + 1), at least 
in the high temperature regime fc^T 3> (i. For low mo- 
menta such that Ek < /x, the semi-classical approxima- 
tion correctly reproduces the leading term (feaT/e^) 2 and 
has an error O(l). The relative error is therefore of sec- 
ond order in T. For high momenta fi <C Ek — ksT, 
the semi-classical approximation reproduces the quan- 
tum term with a relative error 0[(P/j,) 2 ]. After summa- 
tion over all k states, one finds for a two-dimensional 
Bogoliubov gas in the thermodynamic limit that both 
the quantum and the semi-classical values of /„ have the 
form: 



2tt< 2 



1 + ln 



k B T 
2(i 



■0[/3 M ln(/3 M )U, (32) 



where £ is the healing length defined by h 2 /m^ 2 = fi. 
These results are summarized in Fig|3j where the semi- 
classical approximation for /„ is compared to the quan- 
tum value as a function of ksT/ \i. 




k B T/^t 

FIG. 3: For a two-dimensional Bogoliubov gas in the thermo- 
dynamic limit, normal fraction f n as a function of the tem- 
perature fcsT. Solid line: quantum prediction. Dashed line: 
semi-classical prediction. In order to have (within Bogoliubov 
theory) a universal function of fcsT/p, we actually plot the 
product of /„ times n£ 2 , the healing length £ being defined 
by h 2 /m£ 2 = fx. 



The last observable that we investigate is the pair dis- 



Within the Bogoliubov approximation, this can be writ- 
ten for a two-dimensional system in the thermodynamic 
limit as: 



2 (2) (r) 



1+1 

n 



d 2 k 

(2tt)2 



cos(k-r) (a^dk + a k a_ k ) (34) 



where n is the total density. For a given k, the expec- 
tation value between square brackets in (|34[) is equal to 



o R . Its value is given by Eq. (fT5)) for the quantum theory 
and by Eq. (|2TJ)) for the semi-classical theory. 

In Fig[4]we plot the pair distribution g^ 2 \r) as a func- 
tion of r for various values of the temperature. The nar- 
row dip which appears in the result of the quantum cal- 
culation originates from the zero-point fluctuations of the 
Bogoliubov modes, and is therefore absent in the semi- 
classical approximation: in the quantum case, the decay 
of the Fourier transform of (r) — 1 at large k is in fact 
algebraic, whereas it is Gaussian in the semi-classical ap- 
proximation. On the other hand, the semi-classical ap- 
proximation reproduces remarkably well the intermediate 
to long-distance behavior already at temperatures as low 
as ksT = 2/x. 




FIG. 4: (Color online) For a two-dimensional lattice Bogoli- 
ubov gas in the thermodynamic limit, pair distribution g^ 2 ' (r) 
as a function of r for different values of the temperature, 
ksT / (i = 0,2,3,5 (from bottom to top). Solid line: quan- 
tum result. Dashed line: semi- classical approximation. In 
the plot, the product of g' 2 ' — 1 with n£ 2 is actually plotted, 
where n is the density, and £ the healing length such that 
h 2 /(m£ 2 ) = (i. For the Bogoliubov gas, this product is in- 
deed a universal function of hsT/fi and r/£. Here the lattice 
spacing is 0.07£. 



C. In the canonical ensemble 

In the language of [3!| , the semi-classical method dis- 
cussed in the previous sections can be seen as a "simple 
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coherent" scheme from which the noise terms have been 
dropped. This suggests that a similar procedure may be 
applied to the "simple Fock" scheme in order to devise 
a semi-classical method that works in the canonical en- 
semble, i.e. at a fixed number N of particles. 

The building block of this scheme is the Fock state 
defined as usual as: 



\N:1>) 



1 



(35) 



|0) is here the vacuum state and the operator creates 
a particle in the (not necessarily normalized) ip state: 



(36) 



By projecting both sides of (p} onto the subspace with 
exactly N particles, it is easy to see that any iV-body 
density operator can be expanded on a family of Fock 
states as: 



Dil>V[if>]\N :ifi{N :if>\, 



(37) 



IMI=i 



where the distribution V is the Fock state equivalent of 
the Glauber-P distribution, and the integral is taken over 
the unit sphere \\ip\\ = 1. The infinite temperature state 
Pcq{T = 0) = 1 is obtained by simply taking a constant 
value for V[ip]. This corresponds to a random selection 
of the wavefunction ip(r — 0) with a uniform distribu- 
tion on the unit sphere \\ip\\ — 1. At finite temperature, 
the distribution function for an interacting gas is 

unfortunately not necessarily regular and positive; as a 
consequence, no stochastic evolution for -0 exists such 
that the thermal density operator p(P) is obtained as the 
average of dyadics of the form \N : ip){N : On the 
other hand, one can find a stochastic evolution ensuring 
that p((3) is the average of dyadics of the slightly different 
form \N : ipi)(N : ipi and ip2 are here independent 
realizations of the Ito stochastic process (39j 



dtp(r) 



dr 
Y 



ho + go- 



N - 1 



Ml 



90- 



N-lEr' dV\^(r')\ 4 



m 



tf(r)+dB(r), (38) 



starting from the common value ip(r — 0), and the corre- 
lation functions of the noise dB(r) satisfy the condition: 



dB(r) dB(r') = 



' 2dV 



Q r Q r , [«5 r , r ^(r)V(r')] , (39) 



where the projector Q projects orthogonally to the ket 

From this exact reformulation of the full many-body 
problem, it is immediate to obtain a canonical version of 
the semi-classical scheme by simply neglecting the noise 
term dB in (f3"8")) . Intuitively this is expected to constitute 



a good approximation of the quantum model at least in 
the high-temperature case, i.e. for 'times' r short enough 
for the effect of the noise terms to remain small. The 
corresponding semi-classical approximation of the den- 
sity operator for the thermal equilibrium state at tem- 
perature T in the canonical ensemble is therefore 



VtP(Q)\N:iP({3))(N :tP((3)\, (40) 



lh«o)||=i 



where ip(f3) has evolved from its initial value "0(0) during 
a 'time' j3 according to the deterministic part of (I38)) . 



d T tp(r, t) 



h +g - 
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|0(r,r)| 2 



-go 
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0(r,r), (41) 



which closely ressembles an imaginary time Gross- 
Pitaevskii equation. 

This semi-classical Fock scheme can be used as the core 
of a numerical Monte Carlo code to study the proper- 
ties of a A^-body Bose gas at thermal equilibrium. From 
the computational point of view, the only non trivial as- 
pect is how to efficiently perform the sampling of V>(0) on 
the unit sphere. The numerical algorithm that we have 
adopted for this purpose is detailed in the appendix [B] 



III. APPLICATION TO THERMAL VORTICES 
IN THE 2D GAS 

In this second part of the paper, we apply the semi- 
classical technique developed in the first part to the study 
of some among the most significant properties of a ho- 
mogeneous two-dimensional Bose gas at thermal equilib- 
rium in the canonical ensemble. This problem of the 2D 
Bose gas is under active experimental investigation. It is 
known theoretically that the 2D Bose gas exhibits the 
Berezinskii-Kosterlitz-Thouless transition [22|, [H, [24j], 
and this transition was recently observed with cold atoms 
in [20l | . An interesting aspect of the experiments with 
atoms is that they have access to vortices [20j, |2l|, so 
that special attention will be paid here to observables 
such as the density and the pair distribution function 
of thermally activated vortices, for which classical field 
methods [9j and in particular the present semi-classical 
field method, are well suited. Our numerical results will 
then be interpreted in terms of simplified analytical mod- 
els, which allow one to unravel the underlying physics. 

The model Hamiltonian used to describe the system is 
the two-dimensional version of the spatially homogeneous 
lattice model (U) with periodic boundary conditions. The 
value of the coupling constant 170 to be used in the calcu- 
lations depends on the details of the atomic confinement 
along the third dimension: here, we assume a harmonic 
confinement in the z direction, with a harmonic oscilla- 
tor length eiho = y^h/mujz much larger than the three- 
dimensional s-wave scattering length a^. In this limit, 
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one is allowed to neglect the energy- dependence of the ef- 
fective two-dimensional coupling constant g [4(| El| , and 
to simply take [Hj 



,9o 



h 2 2V2na 3D 
m ah 



(42) 



To ensure the two dimensional character of the atomic 
gas, we assume that both the thermal energy ksT and 
the mean field zero-temperature chemical potential gon 
are much smaller than the confinement energy fttu z in 
the z direction. Note that the semi-classical approach 
is limited to the weakly interacting gas regime n£ 2 3> 
1, the healing length £ being defined by ft 2 /m£ 2 = 
ngo- Remarkably, this condition reduces to the density- 
independent one mgo/h 2 <C 1 in two dimensions. 
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A. Numerical results 

1. Normal and non-condensed fractions 

The normal fraction (j29]l describes the response of the 
fluid to a spatial twist of the phase [HI, Hg], while the 
non-condensed fraction is simply the fraction of atoms 
in single-particle states other than the zero-momentum 
plane wave / nc = 1 — Nq/N . These two quantities are 
plotted in Fig[5]as functions of the temperature for three 
different values of the interaction strength go, including 
the ideal gas go = 0. The overall behavior is almost the 
same for all the curves: the dependence on temperature 
is always smooth and, as expected, both the normal and 
the non-condensed fractions tend to 1 (0) in respectively 
the high (low) temperature limit. However, whereas the 
shape of the curve giving the non-condensed fraction is 
not qualitatively modified as go grows, the crossover from 
to 1 of the normal fraction turns out to become some- 
how sharper as the interaction strength is increased [43| . 

It is interesting to compare the results for the ideal gas 
case with a (trivial) calculation performed in the grand 
canonical ensemble: as one can see in Fig[5^, the dashed 
line corresponding to the grand canonical prediction sig- 
nificantly deviates from the numerical simulation results. 
A simple explanation for this can be put forward in terms 
of the finite size of the system, which can indeed lead 
to differences between the two ensembles. In particu- 
lar for a Bose condensed ideal gas, the grand canoni- 
cal ensemble predicts unphysically large fluctuations of 
the number of condensate particles [451 l46l. H?) ; although 
this does not significantly affect the normal and the non- 
condensed fractions plotted here, it will have a dramatic 
impact on other quantities like the density fluctuations 
and the mean vortex density that will be studied in what 
follows. 

In order to fully clarify this issue, an exact canonical 
calculation can be performed by means of the standard 
canonization procedure [48]: the analytical predictions 
for the normal and the non-condensed fractions are plot- 



FIG. 5: (Color online) Normal fraction /„ (black) and non- 
condensed fraction / nc (red) as functions of temperature for 
a two-dimensional Bose gas with N — 1000 particles on a 
square box of size L with periodic boundary conditions, (a) 
Ideal Bose gas. (b) Interacting gas with a coupling constant 
go = O.lh 2 /m. (c) Interacting gas with go = 0.333ft 2 /m. 
Symbols: results of semi-classical simulations on a 64x64 grid 
with 2000 realizations. Solid lines: in (a) exact result from the 
canonization procedure (see text); in (b) and (c), a guide to 
the eye. Dashed lines in (a): the grand canonical predictions. 
The temperature is in units of the degeneracy temperature Td 
such that ksTd = 2ixh 2 n/m. 

ted in Fig|5k and compared to the Monte Carlo ones. 
The agreement is remarkable. 

2. Density fluctuations 

In Fig[S] we plot the temperature dependence of the 
pair distribution function (f3"3")l of the gas evaluated at 
coincident points r = r', i.e. g^(0) [49| . In [ll[ this 
quantity was related in a classical field model to the no- 
tion of a quasi-condensate density in the low temperature 
superfluid regime, tlqc — n y2 — g( 2 ) (0). In the figure, 
the dependence of g^ifi) is shown for three values of the 
interaction strength mgo/h 2 = 0, 0.1, 0.333. In the ideal 
gas case go = 0, the Monte Carlo results are in remark- 
able agreement with the exact canonical results obtained 
from the canonization procedure 50] ; on the other hand, 
at low temperatures, when a significant condensed frac- 
tion is present, the grand canonical prediction g^ (0) = 2 
strongly differs from the canonical results and becomes 
physically incorrect. Concerning the dependence on the 
interaction strength go, our simulations confirm the ex- 
pected trend that an increase of the interaction strength 
go at a fixed value of the non-condensed fraction corre- 
sponds to a strong decrease of the density fluctuations. 

Comparing FigJS] to FigJSl it is immediate to see that 
density fluctuations are already significant in the range of 
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temperatures corresponding to the rapid increase of the 
normal fraction. This shows that density fluctuations 
may play an important role in the superfluid transition 
of a 2D gas [T3, [II]. 
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FIG. 6: (Color online) Pair distribution function <? (0) as a 
function of temperature for the same parameters as in FigO 
Symbols: results of the semi-classical simulations. From top 
to bottom, the value of the coupling constant increases from 
go = (black stars) to go = 0.1ft 2 /m (red diamonds) and 
0.333ft 2 /m (green circles). Solid lines: for go — the exact 
result from the canonization procedure, for go > a guide to 
the eye. Horizontal dashed line: grand canonical prediction 
g^ 2 \0) = 2 for the ideal gas. The temperature is in units of 
the degeneracy temperature T d such that ksT d = 2nh 2 n/m. 



3. Vortex density 

In the semi-classical theory, it is straightforward to de- 
fine a vortex density by looking for the vortices that ap- 
pear in each stochastic realization of the classical field 
ip(r). This is an advantage with respect to e.g. Path 
Integral Quantum Monte Carlo methods [l2T |. 

The field ^(r) of the semi-classical method, initially 
defined on a lattice, may be extended to any point of the 
continuous space by means of the Fourier formula 

^(r) = ±]r ak e* kr , (43) 

k 

where the ak are the Fourier components of the field on 
the lattice. As usual, vortices correspond to nodes in 
the field ip with a non-zero circulation; numerically, they 
can be efficiently and precisely located by calculating the 
circulation of the phase gradient of the field tp around 
plaquettes of much smaller size than the original lattice 

cell mi. 

Numerical results for the mean density of positive 
charge vortices n V) + as a function of temperature for var- 
ious interaction strengths are shown in FigJT^. Thanks 



to the periodic boundary conditions, each realization of 
the field has the same number of positively and nega- 
tively charged vortices, which implies n Uj _ = n Vt +. For 
the considered finite size system, there is no qualitative 
difference between an ideal and an interacting gas: in 
both cases, the vortex density varies roughly linearly with 
temperature at high temperature, while it decreases very 
rapidly at low temperature. Looking at the same data 
on the logarithmic-reciprocal scale of panel (b) , it is easy 
to observe that the low temperature decrease of n Vl + 
roughly follows an activation law of the form oc e ^ A / fc B T . 
A thorough and analytic explanation of this central is- 
sue will be given in section IIII Bl for the non- interacting 
,9o = case and in Sec lIIICI for the general case. 

4- Pair distribution function for vortices 

As a last observable, it is interesting to look at the 
pair distribution function for vortices. In analogy with 
the pair distribution functions for particles in a gas, and 
restricting for simplicity our attention to the case of op- 
posite charge vortices, this may be defined as 

G?i._(r) = (p Ol+ (0K_(r)). (44) 

For a given realization of the field, p v ,±(r) is here the 
sum of Dirac deltas 8{r — y v< ±) centered on the loca- 
tions r Vi ± of the positive (respectively negative) charge 

(2) 

vortices. The angular average of G v ' + _ is plotted as a 
function of the distance r in Fig[H]for different values of 
the coupling constant <?o and temperature. 

In FigEK, a high temperature (but still degenerate) 
case is considered, where both the normal and the non- 
condensed fractions are close to unity: a peak appears 
in all curves at r = as well as a plateau at larger 
vortex separations r. The former is a consequence of 
the effective attraction among opposite charge vortices, 
while the latter corresponds to the decorrelated value 

(2) 

Gy j ~ n v ,+n v ,-- These numerical results indicate a 

weak dependence on the interaction strength, and are in 
good agreement with the known result (not shown) for 
the ideal gas in the grand canonical ensemble [1^, [3(| HH . 

In Fig[5j3, the considered temperatures are low enough 
to be in the regime where n v + drops very rapidly with 
T. For each value of the interaction strength go, the 
temperature is selected to give a roughly fixed vortex 
density. A noticeable difference between the ideal and the 
interacting gas cases appears: the correlations between 
opposite charge vortices have a much longer range in the 
ideal gas than in the interacting one. 

A more intuitive representation of these issues is given 
in FigOU where the locations of the vortices are shown for 
some randomly selected Monte Carlo realizations of the 
field. The high temperature case is considered in (al) 
for the ideal gas and in (a2) for the interacting gas. The 
effect of interactions in the low-temperature regime is 
visible in panels (bl) and (b2): the difference in behavior 
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between the ideal (bl) and the interacting (b2) gas cases 
is apparent, the vortex pairs in the ideal gas being much 
larger. 
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FIG. 7: (Color online) Mean density of positive charge vor- 
tices as a function of temperature for various interaction 
strengths. The parameters have the same values as in Fig[5] 
(a) Linear scale, (b) logarithmic scale for the vortex den- 
sity, reciprocal scale for the temperature. Symbols: re- 
sults of the semi-classical simulation, go = (black stars), 
go = 0.1h 2 /rn (red diamonds), go — 0.333ft 2 /m (green cir- 
cles). Solid lines : the exact canonical result (|46fl for go = 0; 
prediction of the activation law model of Sec lHI Cl for go > 0, 
n V) +/n — Cc" A ( T '/' ! £ ,T i with the prefactor C taken as a con- 
stant and fitted to the data (C = 0.134 for go = O.lh 2 /m and 
C — 0.3355 for go = 0.333ft 2 /m). Dashed line: grand canoni- 
cal result for go = 0. Dot-dashed line: Bogoliubov prediction 
for go = for T/Td < 0.15, essentially indistinguishable from 
the solid line in (a) . Note that the circle with the largest value 
of Td/T corresponds to ksT/ngo ~ 1.4, which is on the limit 
of the validity of both the semi-classical field method and of 
the simple model of section [III CI calculating A. 
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FIG. 8: (Color online) Results of the semi-classical simula- 
tions for the angular average G^ 2 \ (r) of the pair distribu- 
tion function for opposite charge vortices as a function of the 
distance r between the two vortices. The parameters have 
the same values as in Fig[5] (a) High-temperature, non-Bose 
condensed regime, temperature T/Td = 2.5/(27r) ~ 0.398, for 
mgo/h 2 — (black stars), 0.1 (red diamonds), 0.333 (green 
circles). The solid lines are a guide to the eye. Horizon- 
tal dashed lines: square of the mean vortex density n 2 + , 
showing the decorrelation at long distances, (b) Low tem- 
perature, Bose-condensed regime. The temperatures are ad- 
justed to have similar vortex densities for the various val- 
ues of go = (black stars, T/T d = 0.35/(2tt) ~ 0.056, 
leading to n v> + ~ 0.28/L 2 ), go = 0.1h 2 /m (red diamonds, 
T/T d = 0.5/(27r) ~ 0.08, leading to n v ,+ ~ 0.23/L 2 ), 
go = 0.333fi 2 /m (green circles, T/T d = 0.625/(2tt) ~ 0.1, 
leading to — 0.23/L 2 ). The solid lines are a guide to 
the eye. In both panels (a) and (b), the cross at r = gives 

the exact value of G v , for the ideal gas, obtained with the 

canonization procedure. The distance r is in units of L and 
G v j is in units of the squared particle density n . 



B. The effect of Bose condensation on the vortex 
density in an ideal gas: Bogoliubov theory 



To understand the simulation results for the vortex 
density in the non-interacting case, a naive approach is 
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to use the grand canonical ensemble. In this case, the 
Glauber-P distribution for the field is indeed Gaussian, 
so that exact analytical predictions can be obtained for 
the vortex density [H, [3(| : 



(45) 



where = h 2 k 2 /2m, the mean occupation numbers are 
given by the Bose formula, rik = l/{exp[(3(Ek — fi)] — 1}, 
and the chemical potential fi is adjusted to have the same 
density of particles as in the canonical ensemble. 

This prediction is plotted as a dashed line in Fig|7] 
While it is able to correctly reproduce the linear behav- 
ior of the canonical result in the high temperature regime, 
it strongly deviates from it at low temperature: the acti- 
vation law observed in the simulations is then replaced in 
the grand canonical ensemble by a quadratic dependence 
on T. As we shall see in what follows, this deviation is 
due to the presence of a condensate, and is similar to 
the one predicted in [52| for a rotating two-dimensional 
ideal Bose gas in the lowest Landau level. Of course, this 
pathology of the grand canonical ensemble can be elimi- 
nated by a canonization procedure for the vortex density, 
as explained in [52|. We give here only the resulting for- 
mula: 



0v+) c 



in Jo 



J^d9e- i9N 



B(0) 



-i6N 



B{9) 



where the generating function B{9) is written as 
B(0) = l[n k (6) 
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FIG. 9: (Color online) For arbitrary Monte Carlo realizations 
of the field with vortices, locations of the positive charge vor- 
tices (red plus symbols) and negative charge vortices (black 
minus symbols) in the field. Parameters as in some curves 
of FigOD (al) T/T d = 2.5/(2tt) ~ 0.398 for g = 0. (a2) 
T/T d = 2.5/(27r) ~ 0.398 for g = 0.333ft 2 /m. (bl) T/T d = 
0.35/(2tt) ~ 0.056 for g = 0. (b2) T/T d = 0.625/(2tt) ~ 0.1 
for go — 0.333ft 2 jm. Note that the realizations shown in (bl) 
and (b2) are not fully typical since they contain several pairs. 



in terms of a modified Bose law 



n k (9) 



1 



(48) 



As one can see in FigJTl the predictions of this formula, 
are in perfect agreement with the simulation results for 
.90 - 0. 

A physical understanding of the strong suppression of 
vortices in the ideal gas when a condensate is present 
can be obtained by means of the following approximate 
treatment based on the Bogoliubov assumption that the 
fluctuations of the field in the condensate mode are neg- 
ligible. The 2D classical field ip can then be expanded 



4>{r) =ipo + ^2 flk " 



k^O 



L 



(49) 



where the condensate amplitude is fixed to the constant 
value 

^f&y /2 = f y -<;fwy /2 ( 50) 



V L 2 



L 2 
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Here (iVo)Bog is the mean number of condensate par- 
ticles in Bogoliubov theory and the mean number of 
non-condensed particles (SN) Bog in Bogoliubov theory 
is given by 



(MOBog = E 



1 



,PE k 



(51) 



Each of the ak's is a complex random variable with a 
Gaussian distribution [53| : 



-Pk(a) oc e 



(52) 



Since the non-condensed part of the held obeys Gaussian 
statistics, the calculation of the mean vortex density can 
be analytically performed, 



K,+) Bo „ 



E 



E k 



4irH 2 (SN} Bog 



k#0 e f E k-l p -(No) Bos /(SN) I 



(53) 

The prediction of this formula is plotted in Fig]7J as a 
dot-dashed line: the agreement with the exact results is 
good. It is apparent that the dramatic suppression of the 
vortices in the presence of a condensate originates from 
the last factor in Eq. (j53"]) , which is indeed exponentially 
small in the number of condensate particles. One can 
note that a similar factor is involved in the expression for 
the probability to have an empty condensate mode in the 
canonical ensemble. On the other hand, the anomalously 
large vortex density in the grand canonical ensemble can 
be explained by the fact that the most probable value for 
the number of particles in the condensate mode is zero 
in this ensemble. 

Before concluding this section, it is important to re- 
mind that (|53p is an approximate expression. A hrst 
necessary condition for its validity is that a condensate 
is present, which implies N ^> (5N) Bog . For a large box 
L ^> A t h (Ath is here the thermal de Broglie wavelength 
Aj h = 2irh 2 /mksT), this condition corresponds to 



(54) 



nA t \ > 2 log(X/A th ). 



Another necessary condition for the validity of (|53|) is 
that the configurations of the field with vortices are still 
well described by the Bogoliubov model originally de- 
rived for a vortex free field. More precisely, Eq. (|50|) has 
to hold also in presence of vortices, e.g. one has to require 
that the mean number of non-condensed particles condi- 
tioned to the presence of a vortex, say in r = 0, remains 
very close to (SN) Bog . This conditional non-condensed 
number is defined as 



(SNY 



W(r = 0)]E M ok 



W(r - 0)]> 



(55) 



where the expectation value is taken over the exact field 
distribution, S is the two-dimensional Dirac distribution 
and the a^'s are the Fourier components of the field. Cal- 
culating (|55p within Bogoliubov approximation leads to 



the validity condition 



(6N) B ™ d - (SN) Bog = (2 



E 



k^O 



(^o)Bog 

(SN) Bog 
1 



1 x 



of3E k _ 1 



(SN) Bog 



« (SN) Bog . (56) 



In the large box limit L 3> A t h, this condition reduces to 
the simple condition 



4tt 2 

nX 2 th « — [log(i/A th )] 3 , 



(57) 



where the numerical coefficient A — E q ez 2 * 9 4 — 
6.0268. Note that the two conditions (JHJ) and (57]) are 
well compatible in the large box limit L 3> Ath, and de- 
fine a finite validity interval for the Bogoliubov formula 

(ED. 



C. General analytical model for the vortex density 

In this subsection we provide a physical explanation to 
the numerical observation that the vortex density follows 
an approximate activation law at low temperature. This 
is done by developing a simple and physically transparent 
model whose predictions turn out to be in good quanti- 
tative agreement with the semi-classical simulations pre- 
sented in section UlI A[ for both the ideal and the inter- 
acting cases. 

The idea is to look for an approximate field distribution 
of the form 



PsimpicM =e~^S(N 



(58) 



where \\^p\\ 2 = dVJ2 r IVK 1 ")! 2 ; with a suitably chosen en- 
ergy functional U[ip]. As a temperature independent en- 
ergy functional (e.g. the Gross-Pitaevskii one ((Sj)) would 
introduce an unacceptable cut-off dependence [551 ] . we 
are forced to allow for a temperature dependence of U. 

In the ideal gas case, we can reproduce the reasoning 
of Sec lII CI starting from a different representation of the 
infinite temperature density operator, 



p(r = 0) = Dip 



-\W\ 2 
N\ 



\N:i/))(N:^\, 



(59) 



which comes from the projection of the standard over- 
completeness relation for the Glauber coherent states 
onto the iV-particle subspace. Note that ip now runs over 
the whole functional space and is no longer restricted to 
the unit sphere. The evolution |]4"T]) is then applied to 
each initial Fock state; in the go = case, this can be 
solved analytically. Taking the field ip at 'time' (3 rather 
than at time as integration variable, we can write 

pi/3) = I V^P [^\\N : i)/U\\){N : V/IMI I, (60) 
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with the field distribution Po[ip] equal to 



p m 



-\w\ 2 



2 A 



E k |a k | 2 (e^-1) 



(61) 



at is here the Fourier component of the field ip on the 
normalized plane wave 

e ik-r/yi/2_ The ii^n dependent 

prefactor allows for fluctuations of \\ip\\ 2 at most of or- 
der 0(N 1 ^ 2 ) around AT, which, in the large N limit, is 
a relatively small quantity as compared to N. By ap- 
proximatin g th e prefactor with a Dirac delta imposing 
\\ip\\ 2 = N (56|, we finally obtain the desired form ([58]) . 
with the energy functional 



W] = J] |a k | 2 fc B T( ( 



,0E k 



!)• 



(62) 



For the eigenmodes of energy E k <C k B T, this en- 
ergy functional essentially reduces to the non-interacting 
Gross-Pitaevskii energy functional, while for the eigen- 
modes of energy Ek 3> fcsT the large value of e@ Ek 
strongly reduces the modulus of Ok, as required by the 
Bose law for a quantum field. 

This construction can then be heuristically extended 
to the interacting case. Restricting ourselves to relatively 
high temperatures fc^T 3> gon, we can assume that the 
modes for which the interaction energy plays a significant 
role have an energy < gon and can be treated within a 
classical field treatment. This amounts to adding the 
usual interaction term of the Gross-Pitaevskii energy 
functional [57| to the ideal gas functional 



E 

k 



\a k \ 2 k B T(e^ - 1) + | J d 2 v H 4 . (63) 



As the norm of ip is fixed to N in (|55|) . the energy func- 
tional U can be rewritten in the more convenient form 



N 



£K| 2 fcsT(e^-l) 



goN 2 



which is invariant under multiplication of ip by a global 
factor, and allows to formally relax the condition \\ip\\ 2 = 
N. 

The fact that the formation of vortices at low temper- 
ature is an activated process results from the fact that 
the minimal value of U[ip] for a field with at least one 
node is strictly larger than the absolute minimum of U[ip] 
(which corresponds to a nodeless ip). The activation en- 
ergy A(T) is given by: 



A(T) 



min U[tp] — min ^[V 1 ]? 

ij) with a node ip nodeless 



(65) 



and its temperature dependence originates from the tem- 
perature dependence of the energy functional U . In the 
regime k B T -C A(T), the probability to have the field 
with at least one node has the activation form: 



Pnodc 



D -A(T)/k B T J^withanodc r 



ip nodeless 



where the fraction in the right-hand side has an entropic 
origin and is expected to be a slowly varying function of 
T. 

The general strategy to calculate A is what follows. 
Assuming without loss of generality that the node is in 
r = 0, the k = Fourier component ao of the Bose field 
can be expressed in terms of the other components: 



a = 



(67) 



The energy functional U[ip] is then a function of the ak/o 
only and can be minimized without having to impose any 
further constraint. 

The calculation of A(T) is straightforward in the ideal 
gas case. We have to impose that the first order differ- 
ential of U[ip] with respect to the ak's vanishes, which 
leads to the condition 15a 



A/N 



A/N - m ' 



(68) 



where rjk = kBT(e l3Eh — 1). Inserting this equation into 
(fBTf gives a closed equation for the activation energy, 



A/N 



k^O 



Vk-A/N' 



(69) 



A graphical reasoning shows that there exists a unique 
solution in the interval < A/N < r)2v/L, which is the 
smallest root of Eq. (|69|) and thus gives the value of A. 
In the large box limit L 3> Ath, one has the analytic 
expansion: 



A = 



N 



(Emo n k 1 ) 2 



d 2 * IV 7 ! 4 ; whose leading term reduces to 

(64) TTh 2 



mlog(L/A t h) 



(70) 



(71) 



Remarkably, the condition to be in the activation regime 
A ^> k B T is equivalent to the condition (|54|) for Bose 
condensation, N ^> ( (5 TV) Bog- Note also that the lead- 
ing term in (|70[) coincides with the activation part of the 
Bogoliubov result (j53"|) . The successive term gives a cor- 
rection to A which is negligible as compared to k B T pro- 
vided that the validity condition Q57p for the Bogoliubov 
theory is satisfied. 

In the interacting case, a numerical minimization of 
U[ip] in the subspace of the fields with a node in r = 
is performed with the conjugate gradient method. As an 
initial guess, a ip with random complex Fourier coeffi- 
cients ak^o is used. We find that the minimizing field 
ipo has a uniform phase and has a double node in r = 0. 
This means that ipo may be taken real and corresponds 
to the superposition of two, oppositely charged vortices 
located in the origin. 
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Note that it is possible to reduce the energy U by 
continuously transforming this field configuration into a 
nodeless configuration with just a dip in the density at 
r = 0. On the other hand, a continuous transformation 
of this field configuration into a configuration with a pair 
of closely spaced opposite charge vortices corresponds to 
an increases of the energy U. 
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FIG. 10: (Color online) (a) Cut along a;-axis of the field 
tpo minimizing the energy functional U[ip] over the fields 
with a node at the origin. Black solid line (the broadest 
hole): g = 0, T/T d = 0.35/(2tt) ~ 0.056; red solid line: 
go = 0.lh 2 /m, T/T d = 0.5/(2tt) ~ 0.08; blue solid line (the 
narrowest hole): go = 0.333ft 2 /m, T/T d = 0.625/(2tt) ~ 0.1. 
The total number of particles is N = 1000. The dashed lines 
for go > correspond to a field value (m/So) 1 ^ 2 , where fi is 
the Lagrange multiplier defined in Eq. (|74[) . (b) For a semi- 
classical Monte Carlo realization of the field with a single vor- 
tex pair with a small radius, comparison of the density profile 
of the field (green solid line) with the one of the minimizer tpo 
of U[ip] with a node (black solid line). Here go = 0.333ft 2 /m, 
T/T d = 0.5/(2tt) ~ 0.08, the vortex pair diameter is ~ 0.03L 
and the origin of the coordinates was redefined to match the 
location of the vortex pair. 

In FigfTUb we show a cut of the field tpo along x axis 
for the same parameters as in FigJSJa. In FigllOb we com- 
pare the corresponding density profile to the one of a ran- 
domly chosen Monte Carlo realization with a small radius 
vortex pair: there is an acceptable agreement, specially 
considering the significant density fluctuations in the sim- 



ulation result even at the low value of the temperature 
considered here. It is apparent on FigfTUk that the field 
ipo has a slowly varying long-distance tail in the ideal gas 
case, whereas it rapidly reaches its limiting value in the 
interacting case. This can be understood analytically as 
follows. 

For the ideal gas in the thermodynamic limit, one uses 
Eqs. ([57)l and (frj5|) , neglecting A/N with respect to rjk 
(for k > 2tt / L) and then replacing the sum over k by an 
integral, to obtain the approximate expression 



Mr) ^ a ° L J^ 



<i 2 k 1 — cos k • r 



(27T) 



Vk 



(72) 



which holds for r much smaller than the box size L. In 
the limit of large r 3> Ath, the integral is dominated by 
the contribution of the low momenta, which results in 
the functional form 



ip {r) tx ln(r/A th )- 



(73) 



In the interacting case, a sort of generalized Gross- 
Pitaevskii equation can be derived, expressing the fact 
that ipo is an extremum of U[ip] under the constraint 
that the norm is constant and a node is present in r = 0, 



k B T ( e -^ 2 v 2 / 2m _ ij + ffo |^|2 _ M j Mr) 

t iLa + g I |^o| 2 ^o ) S(t). (74) 



/i is here the Lagrange multiplier associated to the con- 
dition of a constant norm for ip. Using the numerical 
fact that tpo is a real function and assuming that at 
large distance from the origin the laplacian term V 2 ?/>o 
is negligible, it is easy to see that nas to converge 
to the limiting value /J,/ go- The normalization condi- 
tion ||V'o|| 2 = N then leads to /i ~ g^n in the large L 
limit. To see how fast ipo reaches its limiting value, we 
set ipo(r) = (^/g ) 1 / 2 [l + i^(r)] and we linearize the equa- 
tion in tp, 



-/3ri 2 V 2 /2m 



1 



2/x 



ip(r) ~ 0. 



(75) 



We heuristically assume that, at large r, tp varies slowly 
at the scale of the thermal de Broglie wavelength. The 
first operator in the above equation may then be approx- 
imated by the usual kinetic energy operator, so that 



ft 2 V 2 
2m 



2[i 



tp(r) ~ 0. 



(76) 



The solution is tp(r) cx Ko(2r/£) where £ is the healing 
length, and K (u) is a Bessel function that tends to zero 
at large u as e~ u /it 1 / 2 . As a consequence, at large r, 



Mr) = ~ 



1/2 



l + O 



(e- 2 ^) 



(77) 
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FIG. 11: (Color online) (a) Activation energy A(T) as a func- 
tion of {gori/hsT) 1 / 2 at a fixed particle density n for increas- 
ing system size L/X t h = 6, 12, 24, 48 (thin solid lines, respec- 
tively black, red, green, blue, from top to bottom; the crosses 
are the actually calculated values and the lines are a guide to 
the eye). The dashed line is the upper bound Eq. (|78p for an 
infinite system size. The thick solid line is the improved upper 
bound discussed around Eq. (|79[l . plotted for (ngo/kBT) 1 ^ 2 > 
0.01. (b) Vortex density as a function of the total particle 
number (for increasing system sizes) for fixed values of the 
density n and the temperature T = 0.5 T d /(27r) ~ 0.08 T d , 
and a coupling constant go = 0.333 h 2 /m. Circles: semi- 
classical simulations. Solid line: prediction of the activation 
law QAAe~ A / kBT where the numerical factor 0.44 was fitted 
to the data. 



Since k B T 3> gon, one indeed finds that, at large r, (p(r) 
varies slowly at the scale of A t h, so that our heuristic 
assumption is a posteriori justified. 

This discussion reveals a key difference for the acti- 
vation energy between the ideal gas and the interacting 
gas in the thermodynamic limit. While in the ideal gas 
case the activation energy tends to zero in the thermo- 
dynamic limit, in the interacting case it has a non-zero 
limit. This point is illustrated in FigfTTk. where we plot 
the activation energy A as a function of (gon/ksT) 1 / 2 
for increasing system sizes at a fixed particle density n. 
Away from the origin go — 0, a nice convergence towards 



a universal curve is obtained, while the dependence of A 
on the system size remains apparent for go = 0. A phys- 
ical interpretation of this fact is that, in the interacting 
case, the minimizer ipo exponentially converges to a lim- 
iting value for r ^> £, whereas in the ideal gas case it is 
logarithmically sensitive to the box size L. 

As a consequence of a non-zero value for the activa- 
tion energy in the thermodynamic limit, we expect that 
the vortex density is an intensive quantity for the inter- 
acting gas. This is confirmed by results of Monte Carlo 
simulations for the vortex density as a function of the 
system size at fixed density and temperature: note on 
Fig lllb how the vortex density is remarkably constant in 
the thermodynamic limit. 

As is apparent in FigfTTk. the convergence of the acti- 
vation energy A to its thermodynamic limit value is not 
uniform in ngo/k B T but becomes slower and slower for 
smaller interaction strength. Analytical results can be 
obtained for an infinite size system, as detailed in the 
appendix [C] One finds an upper bound on the thermo- 
dynamic limit value Aqo of the activation energy, 



< 



2Trh 2 n 1 - 2ng /k B T 
m \n[k B T / (2ng )\ 



(78) 



This explicit upper bound is represented by a dashed line 
in FigfTTk. It shows that tends to zero for vanishing 
interaction strength, which makes a physical link with 
the ideal gas result Eq. (|7ip in the thermodynamic limit 
L/Xth -> oo. 

A better upper bound, though requiring some numer- 
ics, is obtained by performing a variational calculation, 
based on the thermodynamic limit of the ansatz 



1 — cos(k • r) 



exp(E k /k B T eS ) - 1 



(79) 



where J\f is a normalisation factor. The two variational 
parameters are an 'effective' temperature T c q and a > 0. 
The physical motivation for this ansatz, as well as the 
way to implement it in the thermodynamic limit, are 
given in the appendix [Cl The prediction of this ansatz 
is shown as a thick solid line in Fig lllh : it is almost in- 
distinguishable (on the figure) from the numerical results 
for the largest system sizes, except in g = where the 
numerical results suffer from finite size effects. 

The success of this ansatz is due to the fact that it 
reproduces in a fairly accurate way the spatial shape of 
tpo both at short and long distances: In the limit ngo <C 
k B T the energy minimisation leads to T c g ~ T and a ~ 
1.5ngo/k B T. At distances r <C £ one is then allowed to 
neglect a in the denominator of (|79[) . In this way, one 
recovers the ideal gas result ((75)) and, in addition, one 
obtains the normalization factor which depends on the 
interaction strength, 



Mr) 



21n(r/A th ) 
ln(l/a) ' 



(80) 
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In the large r limit r 3> £, the ansatz reproduces the 
exponentially fast convergence of i/'o towards its limiting 
value, Eq. (|77| . with a decay length differing from the 
exact one by a numerical factor close to unity, ~ 1.15. 

From Eq. (|80)) . it is possible to estimate the half-width 
at half maximum of the hole in the density profile ip 2 ,: 

in the g limit, a result growing as A t h i^/^th) 1 ^ 
is found. This prediction is in good agreement with the 
numerical results of FigfTTb for g > and the largest 
sample size, L = 48Ath- 



IV. CONCLUSIONS 

In this paper, we have introduced a semi-classical field 
method for the study of the thermal equilibrium state 
of an ideal or weakly interacting Bose gas at finite tem- 
perature. We have validated the method by verifying 
that it does not suffer from ultraviolet divergences and 
it provides quantitatively accurate predictions as long as 
the temperature is higher than the chemical potential 
of the gas. The method being based on a probability 
distribution in the functional space of c-number wave- 
functions, it appears as being particularly well suited to 
the study of thermal vortices, in contrast to standard 
Quantum Monte Carlo techniques. 

As a first application of the method to a system of 
current experimental interest, we have calculated in this 
paper the density of thermal vortices in a spatially ho- 
mogeneous, two-dimensional Bose gas at thermal equilib- 
rium and we have characterized the spatial correlations 
between the positions of opposite-charged vortices. The 
numerical results are then used as a starting point to de- 
velop simple analytical models and obtain an insight in 
the physics of the system in the different regimes. 

In both the ideal and the interacting cases, in the 
low temperature limit, the vortex density depends on 
temperature according to an activation law of the form 
exp(— A/ksT), with an activation energy A weakly de- 
pendent on temperature. For the ideal gas, A is non-zero 
for a finite size system, because Bose-Einstein condensa- 
tion takes place in such a system at low enough tem- 
perature; for the same reason, A depends on the system 
size and tends logarithmically to zero in the thermody- 
namic limit. For the interacting gas, A has a non-zero 
value in the thermodynamic limit, reached for a system 
size larger than the healing length £; this thermodynamic 
limit value of A tends to zero logarithmically in the limit 
of a vanishing interaction strength. 

Finally, we have studied the spatial correlations be- 
tween the positions of vortices. At high temperatures, 
no qualitative difference appears between the ideal and 
the interacting cases. On the other hand, at low tempera- 
tures (i.e. in the activation regime), the correlations have 
a much longer range in the ideal gas, which corresponds 
to the existence of larger size vortex pairs. 
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APPENDIX A: QUANTUM AND 
SEMI-CLASSICAL GLAUBER-P 
DISTRIBUTIONS FOR A SINGLE BOGOLIUBOV 
MODE 

In this appendix we calculate the exact Glauber-P 
distribution P(j) and its semi-classical approximation 
-Psc(7) for the thermal density operator of a single mode 
Hamiltonian of the form ([16]) . 

The imaginary-time evolution of the Glauber-P dis- 
tribution P(j) is very similar to the one ^ of the full 
many-body Hamiltonian: 

d T P(j) = -E 1 ( 7 )P( 7 ) 

- {d, [Fx ( 7 ) P( 7 )] + f d 2 P( 7 ) + c.c.}. (Al) 

In the (x, y) variables defined as the real and the imagi- 
nary parts of the field 7 = x + iy, the mean-field energy 
£1(7) and the drift force have the simple form: 

£1(7) = (E k + 2/i) x 2 + E k y 2 (A2) 

£1(7) = -\[{E k +2n)x + iE k y], (A3) 



while the diffusion matrix is non-positive definite due to 
the squeezing terms c 2 and (c^) 2 in the Hamiltonian. 

The analysis of the exact P( 7 ) is most easily done by 
looking at its Fourier transform, i.e. the normally ordered 
characteristic function 1331] : 



xp(0 



(A4) 



where the expectation value is taken on the normal- 
ized thermal density operator. For a normalized Gaus- 
sian density operator originating from the imaginary- 
time evolution under a quadratic Hamiltonian such as 
(fl"6]) . Wick theorem implies that: 



Xp(£) = exp 



2 (ctct)+r 2 (cc)-2|e| 2 (ctc» 



(A5) 

From the Gaussian structure of X-p(£)j it is immediate 
to see that the Glauber-P distribution is positive and 
regular if and only if: 



c f c) > 



(A6) 



Applying this to the Bogoliubov theory provides the con- 
dition on the thermal mode occupation: 



1 



of)e k 



1 

T>2 



E k + 2y 
E k 



1/2 



1 



(A7) 
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which is plotted in FigQ] 

We now turn to the semi-classical approximation. The 
solution of the evolution of 7 under the drift force, 
Eq. pTj) . is a simple scaling transformation: 



V(0) 



-E k P/2 



1/(0). 



(A8) 
(A9) 



An explicit form of the weight W(x(0), 2/(0); (3) is then 
obtained by inserting the explicit solution (|A8flA9|) into 
(TT2")) and integrating it. The result is a Gaussian distri- 
bution as a function of (x(Q),y(0)): 

W(f3) = cxp{-[(l - e-^ Ek+2 ^)x(0) 2 

+ (l-e-^)y(0) 2 ]}. (A10) 

The semi-classical approximation (| 13[) to the (unnormal- 
ized) Glauber-P distribution is finally obtained by simply 
writing (|A10|) in terms of the final variables (x(j3), y(/3)) 
for which 7 = x(/3) + iy(f3). As the Jacobian of the rescal- 
ing transformation (|A8IIA9[) is a constant (independent 
of x(0) and y(0)), the result is the Gaussian distribution 
(flTf with the widths given by l|20|2ip . 



APPENDIX B: NUMERICAL ALGORITHM IN 
THE CANONICAL ENSEMBLE 

At r = 0, a wavefunction ip(0) has to be randomly 
selected on the unit sphere, and then let evolve until 
t = (3 = 1/ksT according to (l4Tj) . This provides the 
final value °f the wavefunction to be used in (j40| . 

The observables are then computed as averages over the 
different realizations. As (|4Tj) is purely deterministic, this 
reduces to an averaging over the possible initial wavefunc- 
tions i(j(0)- 

In order to improve the statistical properties of the 
Monte Carlo code, an importance sampling technique [6l| 
has been implemented in terms of an a priori probability 
distribution Q[i/j(0)]. The expectation value of a generic 
operator O is rewritten as: 



(6) 



1 

z 



Xty(O)QMQ)] 



11-011=1 



(N : WMO\N : ^{(3)) 

QMO)} 



( B1 ) 

where Z is the normalisation factor. If the distribution 
of the initial wavefunction "0(0) is sampled with a prob- 
ability law proportional to Q[ip(0)}, one is left with the 
average of a quantity 



(N : tP(P)\6\N : 0(/3)) 



(B2) 



which can be made flatter by means of a clever choice 
of Q[V>(0)]. This provides significant improvement to the 
statistical properties of the Monte Carlo code. In our 
simulations the form 



Q[i>(0)} = (N:i>((3)\N:Tp(P)) = \\ip((3)\ 



2JY 



(B3) 



is used for the a priori probability distribution Q[ijj(0)]. 
This choice was motivated by the requirement that the 
integrand (|B2|) be flat at least for the calculation of the 
partition function, i.e. the trace of the density opera- 
tor. In the numerical code, the sampling of Q[V>(0)] is 
performed by means of a standard Metropolis algorithm 
based on rotations in the single-particle Hilbert space, 
similarly to what was done in [62j . 

Another way of sampling Q[V>(0)] (not used in this 
work) could be the following. At each step of the 
Metropolis algorithm, a multiplication of the amplitude 
of V'(O) on one (randomly chosen) mode of the system by 
a random complex number z — e x e la is proposed. The 
phase a is uniformly distributed in [0, 2k[ and the log- 
arithm A of the modulus has an even probability distri- 
bution over the real axis. Subsequently one renormalizes 
■0(0). One can check that this procedure preserves the 
detailed balance condition required by the Metropolis al- 
gorithm, since the probability distributions of z and of 
1/z coincide. 



APPENDIX C: UPPER BOUND ON THE 
ACTIVATION ENERGY A IN THE 
THERMODYNAMIC LIMIT 

In this appendix we derive an upper bound on the ther- 
modynamic limit value of the activation energy Eq. (|65| 
for an interacting gas go > 0. 

To take the thermodynamic limit in the energy func- 
tional U[ip], we set 



0(r)=A/7(r) 



(CI) 



where /(0) = and f(r) reaches rapidly unity at large 
r/£. The normalization factor is given by 



(l/l 5 



(C2) 



where we have subtracted and added one to \f\ 2 . As 
|/| 2 — 1 is an exponentially narrow function of r for 
L — * 00, the integral in Eq. (|C2|) rapidly converges in 
the thermodynamic limit, so that we get the expansion 



|AA| 2 = n 



1 - 



1 

V- 



(|/| 2 -1) + 0(1/L 4 ) 



(C3) 



where the integral is now over the whole plane. This 
allows to calculate the deviation SU^f] between U[ip] 
and the nodeless ground state energy N 2 g /2L 2 in the 
thermodynamic limit, as a functional of /. For Z7o['0] the 
knowledge of the leading order term of the normalization 
factor Af in (|C3|) is sufficient, whereas for the interaction 
energy the 1/L 2 correction is required. We obtain 



SUootf] = nk B T / / 



-/3R 2 V 2 /2m 



/ 



9on 



(\f\ 2 -iy 



(C4) 
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To easily obtain an upper bound on the thermody- 
namic limit value A m of the activation energy, we re- 
strict to the class C of real and isotropic functions / 
such that < /(r) < 1 for all r. Then (|/| 2 - l) 2 = 
(1-/) 2 (1 + /) 2 <4(1-/) 2 , so that 



Aoo < W[f] = nk B T / / 



1 



+2g a n 2 J (l-/) 2 , V/eC. (C5) 

It remains to minimize the energy functional W[f] over 
the class C, which is conveniently done in the Fourier 
space representation 



J^(k)(l-cosk-r) 

/ i§H k ) 



(C6) 



a writing which ensures that /(0) = and / — > 1 at infin- 
ity for a smooth (real) function u(k). This representation 
leads to 



W[f] = 



(C7) 



Imposing that the functional derivative of this expression 
with respect to u vanishes leads to the choice 



>(k) = 



1 



rjk + 2ng 



(C8) 



One can check, at least for k B T > 2ng Q , that the corre- 
sponding function f m {r) indeed takes values between 
and 1 only, so that it belongs to the class C and it is the 
minimizer of W[f] [59| . This results in the upper bound 

M 



Aoo < W[f m ] = 



2nh 2 n 1 - 2ng /k B T 
m ln[k B T/(2ng )} 



(C9) 



The variational ansatz Eq.(|79p is deduced from 
Eq. (|C8|) by replacing the physical parameters T and 2ng 
by two variational parameters T e ff and ak B T c tf. 
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